#!/bin/csh -f
##############################################################################
#	Script to make postscript graphic output of TISC results with GMT 4.0
#	Daniel Garcia-Castellanos
##############################################################################
#Syntax: 	tisc.gmt.job 'model-root-prj' 
##############################################################################
#Note that:
# -This script uses the bc unix command as a calculator. 
##############################################################################

set prj 	= $1

source $tisc_dir/script/tisc.common.gmt.job


set width = 8


psbasemap -JX$size -R$Region -B$ticklabels\:"x (km)":/$ticklabels\:"y (km)":NseW \
	-Y20.5 -X2.5 -K -P >! $ps
source $tisc_dir/script/tisc.common.topo+drainage.gmt.job
pstext	-JX -R -O -G0 -W255/255/255 -K <<END >> $ps 
	$xtime	$ytime	13 0 1 9	$Timenow Ma
END



#EROSION & DEFLECTION:
if (-r $prj.xyzt || -r $prj.st) then
	psbasemap -JX -R$Region -B$ticklabels\:"x (km)":/$ticklabels\:"y (km)":Nsew -X$horz_shift -O -K >> $ps
if (-r $prj.st) then
echo plotting EROSION RATE...
source $tisc_dir/script/tisc.common.erosrate.gmt.job
endif

if (-r $prj.xyzt) then
echo plotting VERTICAL UPLIFT RATE...
source $tisc_dir/script/tisc.common.upliftrate.gmt.job
endif
if (-r $prj.CMP) psxy $prj.CMP* -JX -R -M -W4/0/0/0 -O -K >> $ps 
pstext	-JX -R$Region -O -K -G0 <<END >> $ps 
#	$xt	$yt	11 0 1 3	erosion/deposition rate 
#(shade) & subs. rate (m/Ma)
END
endif
endif




#PIRACY INDEX:
if (-r $prj.xyzt || -r $prj.st) then
	psbasemap -JX -R$Region -B$ticklabels\:"x (km)":/$ticklabels\:"y (km)":nSeW -X-$horz_shift -Y$vert_shift -O -K >> $ps
if (-r $prj.st) then
echo plotting CHI...
awk '(substr($0,1,1)!=">" && substr($0,1,1)!="#"){print $1,$2,$8;}' $prj.bas | \
	xyz2grd -G$tmp.chi.grd.tmp -I$dx/$dy -R$Region 
grdinfo $tmp.chi.grd.tmp
set chimax = `grdinfo $tmp.chi.grd.tmp | grep z_max | awk '{print $5}'`
set chimax_cpt = `awk -v Lxy=$Lxy 'BEGIN{print Lxy/2}'`

echo Calculating slope along rivers...
set min_slope=.10
set dxh = `awk -v dx=$dx 'BEGIN{print dx/3}'`
set dyh = `awk -v dy=$dy 'BEGIN{print dy/3}'`
awk -v smooth=$min_slope -v dxh=$dxh -v dyh=$dyh 'BEGIN{dxyh=(dxh+dyh)/2}{if (substr($0,1,1)!=">" && substr($0,1,1)!="#"){dist=1e3*sqrt(($7-$1)*($7-$1)+($8-$2)*($8-$2)); if (dist<=dxyh) dist=dxyh; if ($5=="R") print $1,$2,smooth+($6-$9)/dist; else print $1,$2,smooth;}}' $prj.xyw | \
	xyz2grd -G$tmp.slope.grd.tmp -I$dx/$dy -R$Region 
grdclip $tmp.slope.grd.tmp -G$tmp.slope.grd.tmp -Sb$min_slope/$min_slope
grdinfo $tmp.slope.grd.tmp


if (0) then
grdsample $tmp.chi.grd.tmp -G$tmp.chi.grd.tmp -I$dxh/$dyh 
grdinfo $tmp.chi.grd.tmp
grd2xyz $tmp.chi.grd.tmp | \
	blockmean -I$dx/$dy -R$Region -E | \
	awk '{print $1,$2,$6-$5; }' | \
	xyz2grd -G$tmp.chi_contrast.grd.tmp -I$dx/$dy -R$Region
grdinfo $tmp.chi_contrast.grd.tmp

grdsample $tmp.slope.grd.tmp -G$tmp.slope.grd.tmp -I$dxh/$dyh 
grdinfo $tmp.slope.grd.tmp
grd2xyz $tmp.slope.grd.tmp | \
	blockmean -I$dx/$dy -R$Region -E | \
	awk '{print $1,$2,$6-$5; }' | \
	xyz2grd -G$tmp.slope_contrast.grd.tmp -I$dx/$dy -R$Region
grdinfo $tmp.slope_contrast.grd.tmp
endif

echo "Slope being ignored in piracy"
grdmath $tmp.chi.grd.tmp = $tmp.piracy.grd.tmp
grdsample $tmp.piracy.grd.tmp -G$tmp.piracy.grd.tmp -I$dxh/$dyh 
grdclip $tmp.piracy.grd.tmp -G$tmp.piracy.grd.tmp -Sb$dxy/$dxy
grdinfo $tmp.piracy.grd.tmp
grd2xyz $tmp.piracy.grd.tmp | \
	blockmean -I$dx/$dy -R$Region -E | \
	awk '{print $1,$2,($6)-($5); }' | \
	xyz2grd -G$tmp.piracy_contrast.grd.tmp -I$dx/$dy -R$Region
grdinfo $tmp.piracy_contrast.grd.tmp



#grdhisteq $tmp.piracy_contrast.grd.tmp -C10 -D | awk '{if (NR>1) if (x!=$1) print $0; x=$1}' > $tmp.piracy_contrast_levels.tmp 
#awk 'BEGIN{print 0}{print $0; lastlevel=$2}END{print lastlevel}' $tmp.piracy_contrast_levels.tmp > $tmp.piracy_contrast_levels_plus.tmp 
#makecpt -Crainbow -T$tmp.piracy_contrast_levels_plus.tmp -Z > $tmp.piracy_contrast_cpt.tmp
#grd2cpt $tmp.piracy_contrast.grd.tmp -Crainbow -Z -E10 > $tmp.piracy_contrast_cpt.tmp
makecpt -Crainbow -T0/$chimax_cpt/50 -Z > $tmp.piracy_contrast_cpt.tmp
grdimage $tmp.piracy_contrast.grd.tmp -C$tmp.piracy_contrast_cpt.tmp  -JX -R -K -O >> $ps
pstext	-JX -R$Region -O -K -G0 <<END >> $ps 
	$xt	$yt	11 0 1 3	Chi (km)
END

endif

endif
endif



#RIVER PROFILE
echo -n plotting RIVER PROFILE... 

source $tisc_dir/script/tisc.common.river_profile.gmt.job

#psscale -D-.3/1.25/2.1/.2 -L -C$tmp.erosrate_cpt.tmp -O -K -B:"@+sed./eros. (mm/yr)":s >> $ps 


#CHI PROFILE
echo -n plotting CHI PROFILE... 
set zmin = -200
set zmax = 1600 #`echo $zmax | awk '{print $1*1e3}'`
psbasemap -JX-6.8/2.5 -R000/$chimax/$zmin/$zmax \
	-B$ticklabels\:"Chi (km)":/a500f100g1000:"elevation, z (m)":nSE \
	-Y4.5 -X0 -K -O >> $ps

source $tisc_dir/script/tisc.common.chi_profile.gmt.job



#PALETAS DE COLORES:
psscale -C$tmp.topo_palet.tmp 	 -D-1.2/-6.0/14/.2h -B/:"Elevation (m)": -L -O -K >> $ps
psscale -C$tmp.erosrate_cpt.tmp  -D-1.2/-7.1/14/.2h -B/:"Eros.rate (mm yr@+-1@+)": -L -O -K >> $ps
psscale -C$tmp.piracy_contrast_cpt.tmp 	 -D-1.2/-8.2/14/.2h -B/:"Chi contrast (km)": -L -O -K >> $ps

psbasemap -JX -R -Ba0s -O >> $ps


\rm $tmp.*.tmp
